stanEx5
execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}ex12
hierarchical model
class c=1~k
b0c~N(b00,sb0)
b1c~N(b10,sb1)
yc~N(b0c+b1c*x,s)n=50
k=9
x=runif(n,0,20)
c=sample(letters[1:k],n,replace=T)
b00=rnorm(k,10,5)
b0=b00[factor(c)]
b10=rnorm(k,2,1)
b1=b10[factor(c)]
y=rnorm(n,b0+b1*x,5)
qplot(x,y,shape=c,size=I(2))+
scale_shape_manual(values=1:k)
qplot(x,y,col=c)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionestimate as no class
y~N(b00+b10*x,s)ex8-0.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}mdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
fn(mdl,data)## Initial log joint probability = -123310
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22)
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22)
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 28 -161.995 0.000789827 0.00109828 1 1 75
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -162.00
## b0 15.75
## b1 0.83
## s 15.49
## m0[1] 18.81
## m0[2] 20.27
## m0[3] 26.60
## m0[4] 19.26
## m0[5] 18.67
## m0[6] 23.39
## m0[7] 16.82
## m0[8] 27.36
## m0[9] 20.53
## m0[10] 32.04
## m0[11] 15.97
## m0[12] 25.89
## m0[13] 19.34
## m0[14] 28.25
## m0[15] 20.16
## m0[16] 26.65
## m0[17] 29.80
## m0[18] 23.55
## m0[19] 24.53
## m0[20] 22.42
## m0[21] 30.11
## m0[22] 32.20
## m0[23] 23.57
## m0[24] 17.91
## m0[25] 25.88
## m0[26] 17.44
## m0[27] 27.70
## m0[28] 16.13
## m0[29] 17.18
## m0[30] 20.10
## m0[31] 22.97
## m0[32] 23.13
## m0[33] 18.60
## m0[34] 28.39
## m0[35] 23.61
## m0[36] 18.16
## m0[37] 25.54
## m0[38] 17.69
## m0[39] 17.97
## m0[40] 29.07
## m0[41] 23.31
## m0[42] 28.32
## m0[43] 30.32
## m0[44] 30.30
## m0[45] 20.38
## m0[46] 23.12
## m0[47] 24.49
## m0[48] 29.73
## m0[49] 22.21
## m0[50] 19.29
## y0[1] 13.63
## y0[2] 48.27
## y0[3] 7.22
## y0[4] 17.49
## y0[5] 17.68
## y0[6] 36.86
## y0[7] 8.77
## y0[8] 47.78
## y0[9] 22.35
## y0[10] 21.51
## y0[11] 1.30
## y0[12] 26.93
## y0[13] 33.22
## y0[14] 33.51
## y0[15] -6.02
## y0[16] 25.54
## y0[17] 39.50
## y0[18] 27.53
## y0[19] 23.70
## y0[20] 34.32
## y0[21] 47.03
## y0[22] 35.35
## y0[23] 42.43
## y0[24] -4.79
## y0[25] -0.08
## y0[26] 29.59
## y0[27] 18.88
## y0[28] 4.11
## y0[29] 29.39
## y0[30] 26.19
## y0[31] 2.79
## y0[32] 12.46
## y0[33] 33.39
## y0[34] 38.35
## y0[35] 36.89
## y0[36] 1.64
## y0[37] 43.46
## y0[38] 9.08
## y0[39] 30.30
## y0[40] 18.68
## y0[41] 16.52
## y0[42] 4.32
## y0[43] 20.40
## y0[44] 34.81
## y0[45] 32.87
## y0[46] 1.62
## y0[47] 27.84
## y0[48] 34.31
## y0[49] 33.94
## y0[50] 2.92
## m1[1] 15.97
## m1[2] 17.77
## m1[3] 19.58
## m1[4] 21.38
## m1[5] 23.18
## m1[6] 24.98
## m1[7] 26.79
## m1[8] 28.59
## m1[9] 30.39
## m1[10] 32.20
## y1[1] 9.86
## y1[2] 6.20
## y1[3] 39.43
## y1[4] 7.12
## y1[5] 41.92
## y1[6] 1.27
## y1[7] 24.37
## y1[8] 33.12
## y1[9] 28.90
## y1[10] 7.75
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -160.73 -160.44 1.19 0.98 -163.01 -159.44 1.00 718 774
## b0 16.26 16.15 4.31 4.24 9.05 23.14 1.01 312 387
## b1 0.78 0.79 0.40 0.41 0.12 1.43 1.00 414 684
## s 16.25 16.08 1.69 1.71 13.83 19.20 1.00 2313 1543
## m0[1] 19.16 19.14 3.14 3.08 13.83 24.15 1.01 319 433
## m0[2] 20.53 20.56 2.69 2.68 15.95 24.94 1.01 345 505
## m0[3] 26.53 26.52 2.73 2.59 22.00 30.99 1.00 1466 1193
## m0[4] 19.58 19.59 2.99 2.97 14.48 24.35 1.01 324 452
## m0[5] 19.02 19.00 3.19 3.13 13.62 24.08 1.01 318 433
## m0[6] 23.49 23.53 2.24 2.21 19.72 27.09 1.00 652 1119
## m0[7] 17.27 17.21 3.87 3.79 10.75 23.48 1.01 312 387
## m0[8] 27.25 27.25 2.96 2.82 22.46 32.07 1.00 1393 1209
## m0[9] 20.78 20.82 2.62 2.61 16.28 25.08 1.01 354 555
## m0[10] 31.68 31.61 4.75 4.81 24.03 39.48 1.00 794 1233
## m0[11] 16.46 16.35 4.22 4.15 9.42 23.23 1.01 312 384
## m0[12] 25.86 25.86 2.55 2.39 21.52 29.96 1.00 1431 1193
## m0[13] 19.66 19.67 2.96 2.95 14.60 24.41 1.01 325 454
## m0[14] 28.09 28.08 3.25 3.16 22.83 33.35 1.00 1272 1225
## m0[15] 20.43 20.46 2.72 2.70 15.78 24.88 1.01 342 497
## m0[16] 26.57 26.58 2.74 2.60 22.02 31.06 1.00 1463 1193
## m0[17] 29.56 29.51 3.83 3.79 23.34 35.81 1.00 999 1328
## m0[18] 23.64 23.69 2.24 2.19 19.91 27.26 1.00 692 1174
## m0[19] 24.57 24.59 2.31 2.21 20.69 28.34 1.00 1059 1259
## m0[20] 22.57 22.61 2.28 2.26 18.76 26.30 1.01 488 1023
## m0[21] 29.86 29.80 3.96 3.93 23.38 36.34 1.00 957 1290
## m0[22] 31.83 31.76 4.82 4.86 24.10 39.73 1.00 784 1218
## m0[23] 23.66 23.70 2.24 2.18 19.94 27.29 1.00 697 1174
## m0[24] 18.30 18.32 3.46 3.39 12.48 23.83 1.01 314 421
## m0[25] 25.85 25.85 2.55 2.39 21.51 29.95 1.00 1431 1193
## m0[26] 17.85 17.82 3.63 3.56 11.72 23.68 1.01 313 417
## m0[27] 27.57 27.55 3.07 2.96 22.61 32.49 1.00 1350 1241
## m0[28] 16.62 16.51 4.15 4.07 9.70 23.29 1.01 312 381
## m0[29] 17.61 17.58 3.73 3.66 11.30 23.60 1.01 312 406
## m0[30] 20.38 20.41 2.74 2.71 15.71 24.85 1.01 341 494
## m0[31] 23.09 23.15 2.25 2.22 19.31 26.70 1.00 567 1023
## m0[32] 23.24 23.28 2.24 2.22 19.46 26.84 1.00 597 977
## m0[33] 18.95 18.94 3.21 3.14 13.51 24.05 1.01 318 437
## m0[34] 28.22 28.21 3.30 3.20 22.88 33.53 1.00 1249 1225
## m0[35] 23.70 23.73 2.24 2.18 19.98 27.33 1.00 708 1188
## m0[36] 18.54 18.54 3.36 3.27 12.88 23.91 1.01 315 425
## m0[37] 25.52 25.51 2.47 2.32 21.34 29.44 1.00 1438 1281
## m0[38] 18.09 18.08 3.54 3.48 12.14 23.78 1.01 313 421
## m0[39] 18.35 18.37 3.44 3.35 12.59 23.86 1.01 314 418
## m0[40] 28.86 28.83 3.55 3.45 23.16 34.54 1.00 1115 1241
## m0[41] 23.42 23.46 2.24 2.22 19.64 27.02 1.00 636 998
## m0[42] 28.16 28.15 3.28 3.18 22.85 33.46 1.00 1261 1225
## m0[43] 30.05 30.01 4.04 3.99 23.39 36.63 1.00 935 1257
## m0[44] 30.03 29.99 4.03 3.99 23.39 36.60 1.00 937 1257
## m0[45] 20.64 20.67 2.66 2.65 16.10 25.02 1.01 349 527
## m0[46] 23.23 23.27 2.24 2.22 19.45 26.83 1.00 595 1052
## m0[47] 24.53 24.54 2.30 2.21 20.66 28.29 1.00 1038 1259
## m0[48] 29.49 29.45 3.81 3.77 23.32 35.69 1.00 1008 1328
## m0[49] 22.37 22.42 2.31 2.26 18.54 26.19 1.01 463 894
## m0[50] 19.60 19.62 2.98 2.97 14.50 24.36 1.01 325 454
## y0[1] 19.25 19.35 17.04 16.54 -9.56 46.90 1.00 1887 1788
## y0[2] 20.80 21.21 16.42 15.66 -6.19 47.92 1.00 2020 2024
## y0[3] 26.76 26.86 16.29 16.08 0.12 54.35 1.00 1900 1853
## y0[4] 19.25 18.96 17.00 16.36 -8.24 48.38 1.00 1719 1808
## y0[5] 19.08 19.37 16.85 16.54 -9.69 46.30 1.00 1568 1959
## y0[6] 23.38 23.40 16.11 16.06 -3.52 49.31 1.00 1956 2053
## y0[7] 17.86 17.91 16.77 16.48 -10.07 45.99 1.00 1669 1796
## y0[8] 27.67 26.96 16.53 16.14 0.86 56.08 1.00 1791 1682
## y0[9] 19.74 19.67 16.55 16.00 -8.06 47.63 1.00 1902 1971
## y0[10] 31.37 31.41 16.98 16.86 3.11 59.53 1.00 1865 1689
## y0[11] 16.34 15.98 17.13 17.20 -11.82 43.77 1.00 1441 1643
## y0[12] 26.15 25.84 16.93 16.36 -1.19 54.60 1.00 2027 1918
## y0[13] 19.24 19.26 16.55 16.60 -7.95 46.74 1.00 1840 1835
## y0[14] 28.38 28.64 16.36 15.75 2.53 55.89 1.00 1993 1816
## y0[15] 20.36 20.40 16.78 16.67 -7.53 46.79 1.00 1873 1739
## y0[16] 26.86 26.52 16.52 16.46 -0.23 54.49 1.00 1892 1892
## y0[17] 29.82 29.62 16.76 16.15 3.13 57.92 1.00 1954 1821
## y0[18] 24.07 24.25 16.29 16.11 -3.11 51.05 1.00 1945 1969
## y0[19] 24.76 24.34 16.60 16.56 -1.91 52.91 1.00 2021 1777
## y0[20] 22.10 22.33 16.22 16.70 -5.52 48.30 1.00 2014 1953
## y0[21] 30.35 30.15 16.83 16.59 2.88 58.34 1.00 1904 1955
## y0[22] 32.44 31.90 17.31 16.27 3.84 60.81 1.00 1691 1836
## y0[23] 24.16 24.32 16.56 15.90 -3.63 50.73 1.00 1706 1924
## y0[24] 18.46 18.62 16.68 16.26 -9.47 45.80 1.00 1375 1908
## y0[25] 25.92 25.64 16.49 16.38 -1.03 53.84 1.00 1863 2011
## y0[26] 17.71 17.57 16.69 17.02 -9.34 44.45 1.00 1871 1670
## y0[27] 27.14 26.79 16.54 16.33 -0.34 53.41 1.00 1981 1879
## y0[28] 16.79 16.71 17.02 17.02 -11.02 44.82 1.00 1684 1907
## y0[29] 17.85 17.99 16.66 16.43 -9.77 45.43 1.00 1551 1924
## y0[30] 19.80 20.38 16.46 16.79 -7.80 45.80 1.00 1783 1884
## y0[31] 23.29 22.86 15.75 15.51 -1.89 48.90 1.00 1948 1743
## y0[32] 23.40 23.30 16.33 16.39 -3.05 50.38 1.00 2105 1926
## y0[33] 19.20 19.43 16.38 16.87 -8.45 45.16 1.00 2012 1933
## y0[34] 27.97 28.18 16.30 15.74 0.62 53.93 1.00 1899 1893
## y0[35] 23.88 24.19 17.05 16.79 -4.74 50.85 1.00 1778 1837
## y0[36] 18.12 17.91 17.05 16.85 -9.22 45.99 1.00 1637 1837
## y0[37] 25.47 25.27 16.72 16.68 -1.86 52.40 1.00 2109 1863
## y0[38] 18.23 17.82 17.11 17.23 -9.61 45.97 1.00 1572 1563
## y0[39] 18.54 18.56 16.15 16.07 -7.86 44.96 1.00 1968 1931
## y0[40] 28.78 28.25 16.73 15.93 2.16 56.99 1.00 1923 1971
## y0[41] 23.64 23.61 15.97 15.43 -3.08 49.21 1.00 1826 1761
## y0[42] 28.40 28.59 16.71 17.06 1.11 55.35 1.00 2004 1852
## y0[43] 29.94 29.64 17.02 16.68 2.75 58.31 1.00 1953 1865
## y0[44] 29.52 29.71 16.93 16.48 1.69 57.36 1.00 1947 1954
## y0[45] 20.25 20.39 17.09 16.78 -8.34 48.28 1.00 1870 1947
## y0[46] 22.50 22.51 16.51 15.75 -4.76 49.41 1.00 2042 1930
## y0[47] 24.86 25.15 16.93 16.70 -2.41 52.36 1.00 1767 1990
## y0[48] 29.51 29.10 17.02 16.94 2.16 57.76 1.00 1853 1845
## y0[49] 22.32 22.30 17.09 17.16 -6.49 50.66 1.00 1795 1980
## y0[50] 19.13 18.77 16.93 16.45 -7.84 47.99 1.00 1497 1851
## m1[1] 16.46 16.35 4.22 4.15 9.42 23.23 1.01 312 384
## m1[2] 18.17 18.17 3.51 3.44 12.27 23.79 1.01 313 421
## m1[3] 19.88 19.90 2.89 2.87 14.93 24.59 1.01 329 464
## m1[4] 21.58 21.60 2.43 2.41 17.42 25.66 1.01 395 734
## m1[5] 23.29 23.33 2.24 2.21 19.53 26.89 1.00 608 999
## m1[6] 25.00 25.04 2.37 2.26 20.99 28.81 1.00 1271 1321
## m1[7] 26.71 26.72 2.78 2.63 22.08 31.26 1.00 1451 1193
## m1[8] 28.41 28.39 3.38 3.28 22.91 33.84 1.00 1208 1258
## m1[9] 30.12 30.07 4.07 4.04 23.38 36.74 1.00 928 1257
## m1[10] 31.83 31.76 4.82 4.86 24.10 39.73 1.00 784 1218
## y1[1] 16.70 17.03 16.87 16.77 -12.66 44.24 1.00 1433 1929
## y1[2] 18.04 17.95 17.20 16.34 -10.53 46.42 1.00 1989 1919
## y1[3] 19.96 20.21 16.48 15.70 -7.69 47.03 1.00 1813 1873
## y1[4] 21.21 20.85 16.10 15.65 -4.89 47.89 1.00 1969 1854
## y1[5] 23.57 24.16 16.27 16.08 -2.80 49.78 1.00 1687 1587
## y1[6] 25.29 25.04 16.44 15.65 -0.89 52.51 1.00 1958 1895
## y1[7] 27.22 27.20 16.34 15.96 1.17 54.29 1.00 1957 1946
## y1[8] 28.20 28.45 16.68 16.22 2.01 55.57 1.00 2004 1948
## y1[9] 30.78 30.85 16.60 16.24 3.07 57.94 1.00 1951 2054
## y1[10] 32.02 31.93 16.69 16.20 5.25 59.00 1.00 1781 1767
estimate as independent class
all b0l,b1l do not have a distribution
class c=1~k
yc~N(b0c+b1c*x,s) ex12-1.stan
data{
int N;
int K;
vector[N] x;
vector[N] y;
array[N] int c;
}
parameters{
vector[K] b0;
vector[K] b1;
real<lower=0> s;
}
model{
for(i in 1:N)
y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0[c[i]]+b1[c[i]]*x[i];
y0[i]=normal_rng(m0[i],s);
}
}mdl=cmdstan_model('./ex12-1.stan')
data=list(N=n,K=k,x=x,y=y,c=factor(c))
mle=mdl$optimize(data=data)## Initial log joint probability = -881.851
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 -94.2744 0.00374625 0.119494 0.8574 0.8574 111
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 -94.2546 0.00130082 0.00425903 1 1 217
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 200 -94.2546 0.000716441 0.00153958 1 1 218
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -94.25
## b0[1] 11.44
## b0[2] 6.49
## b0[3] 3.28
## b0[4] 10.54
## b0[5] 8.73
## b0[6] 14.00
## b0[7] 20.87
## b0[8] 5.91
## b0[9] 12.52
## b1[1] 0.41
## b1[2] 2.27
## b1[3] 2.18
## b1[4] 1.96
## b1[5] 2.09
## b1[6] 0.23
## b1[7] 4.00
## b1[8] 1.50
## b1[9] -0.66
## s 4.00
## m0[1] 16.48
## m0[2] 42.74
## m0[3] 36.22
## m0[4] 12.53
## m0[5] 10.17
## m0[6] 27.51
## m0[7] 13.07
## m0[8] 3.19
## m0[9] 15.31
## m0[10] 18.48
## m0[11] 7.09
## m0[12] 29.99
## m0[13] 38.25
## m0[14] 2.47
## m0[15] 8.97
## m0[16] 16.99
## m0[17] 43.80
## m0[18] 27.95
## m0[19] 16.41
## m0[20] 17.99
## m0[21] 90.42
## m0[22] -0.70
## m0[23] 6.23
## m0[24] 14.19
## m0[25] 16.78
## m0[26] 29.03
## m0[27] 34.75
## m0[28] 6.60
## m0[29] 12.35
## m0[30] 18.47
## m0[31] 26.36
## m0[32] 19.27
## m0[33] 12.85
## m0[34] 17.47
## m0[35] 16.16
## m0[36] 32.55
## m0[37] 23.64
## m0[38] 12.40
## m0[39] 31.59
## m0[40] 42.07
## m0[41] 6.44
## m0[42] 28.69
## m0[43] 32.30
## m0[44] 45.57
## m0[45] 43.27
## m0[46] 26.75
## m0[47] 16.40
## m0[48] 44.95
## m0[49] 14.65
## m0[50] 37.98
## y0[1] 15.46
## y0[2] 45.32
## y0[3] 35.44
## y0[4] 19.34
## y0[5] 7.07
## y0[6] 30.07
## y0[7] 19.67
## y0[8] 0.66
## y0[9] 21.42
## y0[10] 24.81
## y0[11] 6.86
## y0[12] 38.44
## y0[13] 38.93
## y0[14] 3.97
## y0[15] 0.26
## y0[16] 19.08
## y0[17] 48.03
## y0[18] 23.31
## y0[19] 6.80
## y0[20] 17.07
## y0[21] 90.07
## y0[22] 6.58
## y0[23] 7.32
## y0[24] 16.30
## y0[25] 20.15
## y0[26] 30.78
## y0[27] 30.57
## y0[28] 1.58
## y0[29] 10.41
## y0[30] 13.10
## y0[31] 23.98
## y0[32] 20.03
## y0[33] 14.10
## y0[34] 10.48
## y0[35] 15.15
## y0[36] 29.29
## y0[37] 22.81
## y0[38] 9.05
## y0[39] 33.38
## y0[40] 41.46
## y0[41] 8.26
## y0[42] 35.52
## y0[43] 36.35
## y0[44] 51.51
## y0[45] 50.71
## y0[46] 25.54
## y0[47] 22.70
## y0[48] 45.46
## y0[49] 23.30
## y0[50] 33.72mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 1.3 seconds.
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 1.3 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 1.4 seconds.
## Chain 4 finished in 1.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.3 seconds.
## Total execution time: 1.6 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -104.71 -104.47 3.79 3.81 -111.21 -99.12 1.00 699 1158
## b0[1] 11.40 11.40 6.60 6.67 0.44 21.88 1.00 1499 1037
## b0[2] 6.48 6.57 4.21 4.34 -0.47 13.22 1.01 2554 1445
## b0[3] 3.40 3.23 7.75 7.39 -9.48 15.99 1.00 1207 1207
## b0[4] 10.58 10.49 5.51 5.62 1.63 19.64 1.00 1822 1306
## b0[5] 8.80 8.80 3.67 3.58 2.78 14.97 1.00 3405 1698
## b0[6] 14.10 14.10 5.86 5.55 4.30 23.84 1.00 1572 1327
## b0[7] 20.84 20.87 2.87 2.68 15.90 25.64 1.00 3547 1304
## b0[8] 5.89 6.02 4.49 4.79 -1.65 12.87 1.00 2184 1242
## b0[9] 12.59 12.59 4.52 4.49 5.29 20.14 1.00 2125 1367
## b1[1] 0.41 0.40 1.28 1.25 -1.62 2.52 1.00 1599 1310
## b1[2] 2.28 2.27 0.44 0.43 1.57 3.02 1.00 2449 1728
## b1[3] 2.17 2.18 0.69 0.65 1.04 3.34 1.00 1276 1249
## b1[4] 1.95 1.95 0.41 0.41 1.26 2.62 1.00 1855 1681
## b1[5] 2.07 2.06 0.40 0.40 1.43 2.72 1.00 2655 1881
## b1[6] 0.22 0.23 0.46 0.44 -0.55 0.99 1.00 1570 1310
## b1[7] 3.99 3.99 0.39 0.38 3.37 4.63 1.00 2796 1691
## b1[8] 1.49 1.48 0.39 0.40 0.89 2.11 1.00 2171 1616
## b1[9] -0.67 -0.67 0.37 0.37 -1.26 -0.07 1.00 2019 1596
## s 5.17 5.10 0.65 0.64 4.19 6.33 1.00 1098 1636
## m0[1] 16.48 16.40 2.82 2.76 11.82 21.13 1.00 3148 1573
## m0[2] 42.67 42.68 1.81 1.74 39.68 45.57 1.00 2572 1497
## m0[3] 36.18 36.16 2.65 2.69 31.78 40.58 1.00 2035 1396
## m0[4] 12.62 12.52 5.15 4.99 4.16 21.14 1.00 1234 1281
## m0[5] 10.22 10.22 3.40 3.40 4.85 15.97 1.00 2235 1341
## m0[6] 27.52 27.48 2.07 2.10 24.18 30.91 1.00 2156 1731
## m0[7] 13.10 12.98 5.04 5.15 5.02 21.31 1.00 1824 1365
## m0[8] 3.15 3.20 2.29 2.22 -0.57 6.77 1.00 1953 1173
## m0[9] 15.37 15.37 3.41 3.25 9.52 21.04 1.00 1627 1404
## m0[10] 18.44 18.39 3.98 3.79 11.86 25.06 1.00 1700 1420
## m0[11] 7.08 7.16 4.11 4.24 0.31 13.64 1.00 2553 1446
## m0[12] 30.03 30.01 3.17 3.22 24.84 35.33 1.00 1788 1449
## m0[13] 38.19 38.20 1.89 1.84 34.97 41.21 1.00 2855 1622
## m0[14] 2.43 2.50 2.51 2.46 -1.71 6.45 1.00 1915 1192
## m0[15] 9.01 9.03 2.88 2.86 4.45 13.90 1.00 2207 1442
## m0[16] 17.00 17.02 1.87 1.86 13.96 19.96 1.01 2087 1396
## m0[17] 43.72 43.73 3.38 3.46 38.09 49.15 1.00 2008 1386
## m0[18] 27.96 27.92 2.08 2.12 24.60 31.35 1.00 2145 1735
## m0[19] 16.44 16.42 1.90 1.89 13.31 19.48 1.00 1944 1591
## m0[20] 17.92 17.94 2.28 2.26 14.08 21.62 1.00 2112 1734
## m0[21] 90.25 90.27 4.85 4.81 82.29 98.45 1.00 2268 1606
## m0[22] -0.77 -0.83 3.87 3.78 -7.34 5.53 1.00 1900 1317
## m0[23] 6.23 6.23 2.06 1.97 2.77 9.60 1.00 2103 1648
## m0[24] 14.21 14.16 3.02 2.95 9.18 19.24 1.00 3291 1435
## m0[25] 16.80 16.83 1.79 1.77 13.88 19.72 1.01 2098 1443
## m0[26] 28.98 29.00 2.31 2.19 25.16 32.74 1.00 3379 1375
## m0[27] 34.78 34.76 4.02 4.03 28.29 41.63 1.00 1695 1340
## m0[28] 6.58 6.72 4.34 4.60 -0.71 13.28 1.00 2178 1241
## m0[29] 12.39 12.35 3.22 3.13 7.05 17.74 1.00 3357 1561
## m0[30] 18.46 18.45 2.47 2.45 14.40 22.56 1.00 2427 1520
## m0[31] 26.36 26.36 2.04 2.08 23.03 29.73 1.00 2177 1641
## m0[32] 19.20 19.21 2.18 2.19 15.54 22.70 1.00 2116 1865
## m0[33] 12.79 12.85 3.46 3.47 7.17 18.54 1.00 1754 1325
## m0[34] 17.46 17.48 2.34 2.28 13.55 21.25 1.01 1914 1364
## m0[35] 16.19 16.23 2.14 2.09 12.61 19.71 1.00 1811 1567
## m0[36] 32.50 32.53 2.11 2.06 28.96 35.82 1.00 3208 1392
## m0[37] 23.55 23.46 2.19 2.13 19.94 27.22 1.00 2132 1621
## m0[38] 12.35 12.43 4.23 4.23 5.19 19.20 1.00 1581 1200
## m0[39] 31.54 31.57 2.16 2.10 27.95 34.95 1.00 3248 1372
## m0[40] 42.00 42.06 3.16 3.16 36.67 47.11 1.00 2019 1405
## m0[41] 6.44 6.41 2.09 2.02 2.92 9.93 1.00 2114 1583
## m0[42] 28.58 28.52 2.84 2.78 24.15 33.26 1.00 2193 1926
## m0[43] 32.18 32.15 3.53 3.52 26.62 37.91 1.00 2249 1720
## m0[44] 45.28 45.38 5.09 5.16 37.02 53.47 1.00 2188 1447
## m0[45] 43.19 43.20 1.81 1.74 40.19 46.10 1.00 2537 1474
## m0[46] 26.76 26.73 2.05 2.08 23.42 30.10 1.00 2168 1707
## m0[47] 16.42 16.42 1.91 1.89 13.28 19.48 1.00 1938 1647
## m0[48] 44.98 44.99 4.21 4.03 38.44 52.02 1.00 2223 1653
## m0[49] 14.56 14.56 5.25 5.00 6.09 23.04 1.00 2090 1407
## m0[50] 37.92 37.92 1.89 1.84 34.70 40.94 1.00 2874 1623
## y0[1] 16.60 16.56 5.90 5.52 7.18 26.59 1.00 2227 1974
## y0[2] 42.68 42.60 5.59 5.61 33.52 51.99 1.00 1914 1894
## y0[3] 36.23 36.24 5.88 5.49 26.42 46.09 1.00 2186 1958
## y0[4] 12.38 12.38 7.34 7.31 0.24 24.32 1.00 1438 1519
## y0[5] 10.16 10.16 6.16 6.15 0.06 19.88 1.00 1955 1840
## y0[6] 27.48 27.67 5.66 5.64 18.07 36.59 1.00 2003 1857
## y0[7] 13.01 13.02 7.25 7.08 1.26 25.04 1.00 1948 1918
## y0[8] 3.20 3.12 5.76 5.87 -6.49 12.73 1.00 2035 1768
## y0[9] 15.16 15.20 6.26 6.26 5.23 25.35 1.00 1888 1886
## y0[10] 18.39 18.32 6.55 6.59 7.49 29.71 1.00 1753 1776
## y0[11] 6.99 7.14 6.45 6.46 -3.80 17.46 1.00 1834 1742
## y0[12] 30.11 30.04 6.09 6.07 20.32 40.03 1.00 2107 1895
## y0[13] 38.20 38.17 5.59 5.48 29.11 47.37 1.00 2149 1994
## y0[14] 2.29 2.34 5.86 5.76 -7.47 11.98 1.00 1967 1881
## y0[15] 8.92 8.81 5.79 5.46 -0.49 18.34 1.00 2062 1658
## y0[16] 16.96 16.95 5.54 5.40 7.75 25.97 1.00 1879 1547
## y0[17] 43.86 43.76 6.27 6.02 33.59 54.19 1.00 1862 2013
## y0[18] 27.98 27.93 5.55 5.40 19.10 37.11 1.00 2018 1959
## y0[19] 16.42 16.24 5.47 5.39 7.55 25.62 1.00 2114 2012
## y0[20] 17.94 17.87 5.48 5.29 8.94 27.29 1.00 1964 1732
## y0[21] 90.23 90.05 7.03 6.95 78.71 101.87 1.00 2317 2153
## y0[22] -0.71 -0.65 6.46 6.37 -11.24 10.11 1.00 2042 1878
## y0[23] 6.23 6.20 5.62 5.65 -3.18 15.36 1.00 2077 1440
## y0[24] 14.11 14.04 6.19 6.04 3.79 24.11 1.00 2065 1996
## y0[25] 16.77 16.72 5.47 5.46 7.67 25.57 1.00 1780 1825
## y0[26] 28.88 28.87 5.76 5.68 19.37 38.40 1.00 2047 1962
## y0[27] 34.66 34.63 6.48 6.35 23.86 45.49 1.00 1884 1853
## y0[28] 6.71 6.75 6.64 6.41 -4.37 17.69 1.00 2125 1972
## y0[29] 12.30 12.25 6.20 6.22 2.23 22.14 1.00 2118 1949
## y0[30] 18.51 18.54 5.80 5.79 9.20 28.21 1.00 2089 1881
## y0[31] 26.24 26.13 5.86 6.04 16.58 35.83 1.00 1997 1805
## y0[32] 19.24 19.29 5.66 5.32 9.96 28.60 1.00 2002 1970
## y0[33] 12.80 12.84 6.29 6.32 2.77 23.03 1.00 2046 1986
## y0[34] 17.42 17.36 5.58 5.47 8.34 26.70 1.00 2097 1913
## y0[35] 16.23 16.16 5.80 5.80 6.76 25.73 1.00 1963 1933
## y0[36] 32.47 32.43 5.36 5.29 23.85 41.55 1.00 2091 2042
## y0[37] 23.53 23.54 5.67 5.46 13.89 32.71 1.00 1880 1855
## y0[38] 12.25 12.35 6.65 6.69 1.14 23.07 1.00 1709 1741
## y0[39] 31.54 31.60 5.63 5.72 22.22 40.51 1.00 2305 2012
## y0[40] 42.09 42.33 6.15 6.15 31.94 52.08 1.00 2033 1924
## y0[41] 6.12 6.19 5.60 5.44 -3.16 15.26 1.00 2161 2015
## y0[42] 28.60 28.47 5.98 5.94 19.15 38.09 1.01 1883 1881
## y0[43] 32.09 32.02 6.13 6.10 22.13 42.24 1.00 2303 2001
## y0[44] 45.27 45.39 7.38 7.23 32.60 57.25 1.00 1996 1704
## y0[45] 43.15 43.13 5.41 5.37 34.54 51.99 1.00 1973 1973
## y0[46] 26.64 26.67 5.56 5.47 17.50 35.82 1.00 2013 1893
## y0[47] 16.34 16.38 5.54 5.49 7.40 25.44 1.00 1887 1819
## y0[48] 44.94 44.83 6.61 6.42 33.83 56.10 1.00 2295 1912
## y0[49] 14.47 14.25 7.72 7.49 2.06 27.44 1.00 1981 1416
## y0[50] 37.80 37.74 5.47 5.23 28.55 46.76 1.00 2121 2151
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,shape=c,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
#geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
#geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
#geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
#geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
qplot(x,y,facets=~c,size=I(2),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
estimate as class have relation each other
all b0l,b1l have a distribution
class c=1~k
b0c~N(b00,sb0)
b1c~N(b10,sb1)
yc~N(b0c+b1c*x,s)ex12-2.stan
data{
int N;
int K;
vector[N] x;
vector[N] y;
array[N] int c;
}
parameters{
real b00;
real<lower=0,upper=100> sb0;
vector[K] b0;
real b10;
real<lower=0,upper=100> sb1;
vector[K] b1;
real<lower=0,upper=100> s;
}
model{
b0~normal(b00,sb0);
b1~normal(b10,sb1);
for(i in 1:N)
y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0[c[i]]+b1[c[i]]*x[i];
y0[i]=normal_rng(m0[i],s);
}
}mdl=cmdstan_model('./ex12-2.stan')
data=list(N=n,K=k,x=x,y=y,c=factor(c))
mle=mdl$optimize(data=data)## Initial log joint probability = -266.657
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 -39.1515 0.0339647 538871 0.2692 0.2692 150
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 122 -33.8812 0.000129671 547.855 1 1 178
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.try(print(mle))## variable estimate
## lp__ -33.88
## b00 1.17
## sb0 0.00
## b0[1] 1.17
## b0[2] 1.17
## b0[3] 1.17
## b0[4] 1.17
## b0[5] 1.17
## b0[6] 1.17
## b0[7] 1.17
## b0[8] 1.17
## b0[9] 1.17
## b10 2.50
## sb1 0.45
## b1[1] 4.14
## b1[2] 2.94
## b1[3] 2.60
## b1[4] 2.52
## b1[5] 2.09
## b1[6] 1.57
## b1[7] 5.53
## b1[8] 1.70
## b1[9] 0.29
## s 11.44
## m0[1] 8.91
## m0[2] 31.38
## m0[3] 34.27
## m0[4] 12.22
## m0[5] 2.18
## m0[6] 28.36
## m0[7] 4.43
## m0[8] 5.21
## m0[9] 10.24
## m0[10] 32.07
## m0[11] 1.95
## m0[12] 33.08
## m0[13] 25.19
## m0[14] 5.51
## m0[15] 2.70
## m0[16] 21.84
## m0[17] 44.03
## m0[18] 28.93
## m0[19] 17.83
## m0[20] 14.92
## m0[21] 97.23
## m0[22] 6.89
## m0[23] 3.89
## m0[24] 6.62
## m0[25] 20.40
## m0[26] 12.44
## m0[27] 38.76
## m0[28] 1.96
## m0[29] 4.79
## m0[30] 16.66
## m0[31] 26.86
## m0[32] 16.38
## m0[33] 15.44
## m0[34] 25.15
## m0[35] 16.09
## m0[36] 17.31
## m0[37] 21.35
## m0[38] 10.88
## m0[39] 15.98
## m0[40] 41.80
## m0[41] 3.80
## m0[42] 27.09
## m0[43] 31.21
## m0[44] 37.92
## m0[45] 32.11
## m0[46] 27.38
## m0[47] 17.75
## m0[48] 50.92
## m0[49] 33.55
## m0[50] 24.81
## y0[1] 16.40
## y0[2] 22.61
## y0[3] 43.99
## y0[4] -4.37
## y0[5] 5.40
## y0[6] 11.56
## y0[7] 9.39
## y0[8] -7.71
## y0[9] 11.41
## y0[10] 51.98
## y0[11] 3.20
## y0[12] 11.05
## y0[13] 51.59
## y0[14] 13.30
## y0[15] -11.74
## y0[16] 7.34
## y0[17] 21.90
## y0[18] 45.90
## y0[19] -6.45
## y0[20] 4.72
## y0[21] 108.35
## y0[22] 15.75
## y0[23] -20.95
## y0[24] -0.88
## y0[25] 32.32
## y0[26] 1.66
## y0[27] 54.45
## y0[28] 8.24
## y0[29] 0.68
## y0[30] 32.66
## y0[31] 47.53
## y0[32] 14.66
## y0[33] 20.98
## y0[34] 24.57
## y0[35] 19.70
## y0[36] 21.27
## y0[37] 10.48
## y0[38] 3.96
## y0[39] 28.88
## y0[40] 31.95
## y0[41] -5.74
## y0[42] 37.84
## y0[43] 27.02
## y0[44] 53.68
## y0[45] 28.46
## y0[46] 30.82
## y0[47] -2.94
## y0[48] 49.67
## y0[49] 28.61
## y0[50] 39.75mcmc=goMCMC(mdl,data,wrm=1000,smp=2000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 1 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 1 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 2 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 2 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 3 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 3 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 4 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 4 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 4 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 1 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 3 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 2 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 1 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 3 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 4 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 1 finished in 1.2 seconds.
## Chain 3 finished in 1.2 seconds.
## Chain 4 finished in 1.1 seconds.
## Chain 2 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 2 finished in 1.3 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.2 seconds.
## Total execution time: 1.4 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -128.16 -127.75 4.45 4.36 -136.09 -121.56 1.00 1765 2672
## b00 10.85 10.93 2.65 2.45 6.45 14.92 1.00 4954 5420
## sb0 5.53 5.12 2.54 2.11 2.13 10.16 1.00 2183 1828
## b0[1] 10.22 10.29 3.97 3.77 3.59 16.58 1.00 6858 5653
## b0[2] 8.52 8.60 3.29 3.25 3.05 13.78 1.00 5341 5107
## b0[3] 8.52 8.76 4.59 4.37 0.70 15.47 1.00 5075 4803
## b0[4] 10.88 10.96 3.82 3.60 4.49 17.11 1.00 6748 5749
## b0[5] 9.75 9.76 2.99 2.91 4.83 14.58 1.00 6680 6133
## b0[6] 11.69 11.71 4.12 3.87 5.02 18.54 1.00 5907 5175
## b0[7] 18.93 18.98 3.05 3.03 13.80 23.84 1.00 3737 3172
## b0[8] 8.09 8.25 3.65 3.53 1.93 13.72 1.00 4230 5431
## b0[9] 11.11 11.18 3.46 3.34 5.29 16.68 1.00 6632 5513
## b10 1.54 1.54 0.61 0.54 0.56 2.51 1.00 6480 4133
## sb1 1.66 1.54 0.58 0.46 0.99 2.71 1.00 6736 4916
## b1[1] 0.74 0.73 0.82 0.81 -0.61 2.07 1.00 6883 5807
## b1[2] 2.07 2.07 0.36 0.35 1.49 2.68 1.00 5888 5471
## b1[3] 1.74 1.73 0.45 0.44 1.04 2.48 1.00 5833 4721
## b1[4] 1.93 1.93 0.31 0.30 1.42 2.44 1.00 7114 6069
## b1[5] 2.00 2.00 0.35 0.35 1.41 2.58 1.00 7688 5525
## b1[6] 0.41 0.41 0.34 0.33 -0.14 0.97 1.00 6121 5291
## b1[7] 4.13 4.13 0.41 0.40 3.48 4.82 1.00 4370 3119
## b1[8] 1.33 1.33 0.32 0.32 0.82 1.87 1.00 5195 5527
## b1[9] -0.53 -0.54 0.30 0.30 -1.02 -0.03 1.00 7064 5650
## s 5.06 5.00 0.62 0.61 4.15 6.17 1.00 5593 5711
## m0[1] 17.16 17.16 2.40 2.35 13.24 21.03 1.00 7478 6456
## m0[2] 41.54 41.57 1.85 1.83 38.40 44.52 1.00 6204 5545
## m0[3] 36.17 36.16 2.51 2.44 32.10 40.32 1.00 9774 5518
## m0[4] 15.93 16.06 3.22 3.12 10.43 21.00 1.00 5051 5155
## m0[5] 9.22 9.27 2.62 2.54 4.83 13.46 1.00 6268 6157
## m0[6] 27.71 27.68 1.92 1.88 24.57 30.90 1.00 7107 5251
## m0[7] 13.37 13.44 3.50 3.31 7.55 19.07 1.00 6807 5928
## m0[8] 3.61 3.61 2.25 2.20 -0.07 7.34 1.00 8597 6083
## m0[9] 14.09 14.06 2.48 2.34 10.05 18.24 1.00 6139 5955
## m0[10] 19.86 19.86 3.38 3.33 14.30 25.40 1.00 7080 5864
## m0[11] 9.07 9.14 3.22 3.18 3.72 14.21 1.00 5341 5111
## m0[12] 29.91 29.88 3.11 3.09 24.86 35.08 1.00 7795 5826
## m0[13] 36.90 36.94 1.95 1.93 33.58 40.05 1.00 4972 4535
## m0[14] 3.04 3.03 2.45 2.43 -0.97 7.11 1.00 8506 6209
## m0[15] 8.26 8.30 2.27 2.22 4.42 11.91 1.00 6388 6281
## m0[16] 17.16 17.16 1.86 1.86 14.11 20.19 1.00 8137 5637
## m0[17] 43.63 43.64 3.20 3.19 38.36 48.90 1.00 9451 5870
## m0[18] 28.11 28.09 1.95 1.91 24.92 31.34 1.00 7157 5139
## m0[19] 16.10 16.09 1.73 1.70 13.23 18.89 1.00 7500 5043
## m0[20] 18.84 18.85 2.03 2.00 15.45 22.17 1.00 6902 5863
## m0[21] 90.81 90.82 5.09 4.95 82.39 99.00 1.00 5772 5411
## m0[22] 0.50 0.44 3.56 3.53 -5.27 6.32 1.00 8047 6258
## m0[23] 6.06 6.06 1.84 1.80 3.02 9.14 1.00 8356 5358
## m0[24] 14.97 14.99 2.52 2.46 10.85 19.04 1.00 7126 6251
## m0[25] 16.77 16.78 1.77 1.77 13.88 19.64 1.00 8012 5554
## m0[26] 27.37 27.44 2.44 2.39 23.21 31.26 1.00 3957 3303
## m0[27] 33.72 33.70 3.73 3.64 27.56 39.89 1.00 7685 5999
## m0[28] 8.71 8.86 3.52 3.42 2.73 14.14 1.00 4262 5560
## m0[29] 13.22 13.23 2.65 2.60 8.87 17.47 1.00 6931 6220
## m0[30] 19.45 19.46 2.02 1.99 16.09 22.71 1.00 5750 6319
## m0[31] 26.65 26.62 1.88 1.83 23.57 29.74 1.00 6951 5489
## m0[32] 19.98 19.99 1.99 1.96 16.73 23.26 1.00 7441 5842
## m0[33] 12.76 12.76 2.72 2.67 8.20 17.19 1.00 8731 6005
## m0[34] 18.03 18.03 2.23 2.22 14.39 21.65 1.00 7863 6077
## m0[35] 15.63 15.63 1.80 1.76 12.73 18.58 1.00 7042 5930
## m0[36] 31.01 31.07 2.22 2.18 27.23 34.53 1.00 4200 3636
## m0[37] 23.86 23.87 2.12 2.10 20.44 27.42 1.00 8761 5481
## m0[38] 11.94 11.98 2.89 2.77 7.04 16.59 1.00 7887 6598
## m0[39] 30.02 30.08 2.27 2.24 26.15 33.63 1.00 4123 3702
## m0[40] 41.93 41.93 3.01 3.02 36.97 46.90 1.00 9575 5830
## m0[41] 6.23 6.24 1.85 1.80 3.18 9.30 1.00 8288 5600
## m0[42] 28.35 28.35 2.71 2.69 23.86 32.82 1.00 8533 6366
## m0[43] 31.56 31.57 3.29 3.29 26.15 36.96 1.00 8047 6589
## m0[44] 44.95 44.93 4.94 4.80 36.86 53.08 1.00 9177 6502
## m0[45] 42.08 42.12 1.85 1.82 38.95 45.07 1.00 6381 5525
## m0[46] 27.01 26.98 1.89 1.85 23.90 30.14 1.00 7009 5374
## m0[47] 16.07 16.07 1.73 1.70 13.21 18.87 1.00 7472 5019
## m0[48] 43.62 43.63 3.81 3.77 37.39 49.94 1.00 6826 5018
## m0[49] 15.98 15.98 4.44 4.41 8.75 23.35 1.00 8796 6041
## m0[50] 36.62 36.65 1.96 1.93 33.28 39.77 1.00 4921 4562
## y0[1] 17.17 17.20 5.69 5.64 7.71 26.36 1.00 7950 7366
## y0[2] 41.54 41.58 5.44 5.33 32.52 50.25 1.00 8507 7925
## y0[3] 36.19 36.17 5.70 5.58 26.86 45.53 1.00 8484 7649
## y0[4] 16.00 16.06 6.08 5.91 5.79 25.84 1.00 6841 7117
## y0[5] 9.14 9.13 5.80 5.75 -0.20 18.75 1.00 7931 7594
## y0[6] 27.73 27.72 5.40 5.22 18.78 36.62 1.00 8138 7539
## y0[7] 13.37 13.40 6.13 6.01 3.27 23.41 1.00 7638 7251
## y0[8] 3.56 3.54 5.62 5.47 -5.51 12.83 1.00 8030 7906
## y0[9] 14.21 14.16 5.63 5.59 4.84 23.49 1.00 7756 7567
## y0[10] 19.89 19.88 6.11 6.03 9.80 29.76 1.00 7775 6792
## y0[11] 9.11 9.18 6.03 6.00 -0.76 19.03 1.00 6688 7255
## y0[12] 29.86 29.92 5.95 5.82 20.10 39.66 1.00 7890 7592
## y0[13] 36.97 37.06 5.38 5.31 28.03 45.78 1.00 7755 7337
## y0[14] 3.06 3.03 5.66 5.51 -6.12 12.39 1.00 8069 7334
## y0[15] 8.32 8.29 5.51 5.46 -0.70 17.44 1.00 7733 7628
## y0[16] 17.15 17.08 5.43 5.43 8.37 26.01 1.00 7932 7557
## y0[17] 43.62 43.57 6.06 5.92 33.62 53.56 1.00 8761 7604
## y0[18] 28.06 28.03 5.54 5.36 19.29 37.33 1.00 8136 7810
## y0[19] 16.06 16.08 5.29 5.20 7.38 24.78 1.00 8001 7505
## y0[20] 18.81 18.82 5.46 5.38 9.87 27.78 1.00 7936 7910
## y0[21] 90.76 90.76 7.20 6.97 79.01 102.75 1.00 7013 7182
## y0[22] 0.53 0.56 6.19 6.07 -9.68 10.73 1.00 8004 7789
## y0[23] 6.03 5.99 5.37 5.27 -2.78 14.97 1.00 7943 7643
## y0[24] 14.99 14.92 5.64 5.66 5.78 24.21 1.00 8510 7515
## y0[25] 16.80 16.74 5.36 5.21 8.02 25.68 1.00 8028 7809
## y0[26] 27.44 27.40 5.70 5.56 18.01 36.96 1.00 7021 6967
## y0[27] 33.66 33.69 6.34 6.24 23.39 44.21 1.00 8094 7829
## y0[28] 8.59 8.60 6.24 6.13 -1.66 18.74 1.00 7124 6994
## y0[29] 13.17 13.19 5.67 5.62 3.83 22.31 1.00 7535 7267
## y0[30] 19.43 19.39 5.45 5.35 10.70 28.37 1.00 7289 7584
## y0[31] 26.67 26.68 5.41 5.31 17.83 35.57 1.00 8384 7017
## y0[32] 19.97 19.96 5.48 5.36 11.02 28.97 1.00 8102 7887
## y0[33] 12.74 12.76 5.81 5.71 3.34 22.35 1.00 8199 7693
## y0[34] 17.98 18.00 5.61 5.54 8.81 27.38 1.00 8084 7740
## y0[35] 15.68 15.69 5.42 5.32 6.76 24.49 1.00 8160 7920
## y0[36] 31.04 31.10 5.58 5.48 21.91 40.22 1.00 7240 7178
## y0[37] 23.86 23.84 5.42 5.36 14.87 32.70 1.00 7803 7716
## y0[38] 11.97 12.01 5.92 5.82 2.30 21.64 1.00 8236 8053
## y0[39] 30.03 30.07 5.55 5.41 20.95 39.16 1.00 7440 7646
## y0[40] 42.02 42.10 5.92 5.75 32.26 51.83 1.00 8188 7657
## y0[41] 6.30 6.29 5.49 5.40 -2.80 15.47 1.00 7765 7888
## y0[42] 28.35 28.36 5.74 5.68 19.09 37.60 1.00 8085 8013
## y0[43] 31.64 31.70 5.99 5.93 21.51 41.34 1.00 7802 7412
## y0[44] 44.96 44.86 7.18 7.15 33.18 56.70 1.00 8643 6990
## y0[45] 41.99 41.95 5.45 5.32 33.06 51.00 1.00 7661 7369
## y0[46] 27.05 27.03 5.45 5.38 18.13 35.84 1.00 7599 7636
## y0[47] 16.03 16.07 5.30 5.12 7.27 24.59 1.00 7817 7339
## y0[48] 43.54 43.58 6.37 6.39 32.91 53.83 1.00 8179 7815
## y0[49] 16.11 16.12 6.77 6.68 4.87 27.17 1.00 8543 7976
## y0[50] 36.55 36.65 5.48 5.33 27.44 45.47 1.00 6647 7441
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,shape=c,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
#geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
#geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
#geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
#geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
qplot(x,y,facets=~c,size=I(2),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
ex13
generalized linear mixed model
(X,y)i=1-n
b[b0,b1,...]
linear model y~N(Xb,s)
generalized linear model y~dist.(m=link(Xb),s)
fixed effect b0, b1,...
individual random effect b0+r0i r0~N(0,sr0), b1+r1i r1~N(0,sr1),...
class c=1-k
class effect b0+r0c r0~N(0,sr0), b1+r1c r1~N(0,sr1),...
for y=dist.(m,s)
random intercept model m=b0+r0i+b1*x, m=b0+r0c+b1*x
random coefficient model m=b0+(b1+r1i)*x, m=b0+(b1+r1c)*x
mixed model m=b0+r0i+(b1+r1i)*x, m=b0+r0c+(b1+r1c)*x
note
@ yi=b0+b1*xi+r0i is not useful, ri is included in s
@ yi=b0+b1*xi+r0c is useful, repeated observation can be treated by class effect
@ when appling Poisson dist.(y is larger than 0, error is larger at large x),
but variance is larger than mean (over dispersion),
random intercept model is usefuln=20
x=runif(n,0,10)
r0=rnorm(n,0,1)
y=rpois(n,exp(2+0.1*x+r0))
qplot(x,y)
ex13-0.stan
data{
int N;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
}
model{
y~poisson_log(b0+b1*x);
}
generated quantities{
vector[N] m0;
array[N] int y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=poisson_log_rng(m0[i]);
}
}mdl=cmdstan_model('./ex13-0.stan')
data=list(N=n,x=x,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -3304.39
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 658.086 0.000343231 0.0028631 1 1 17
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 658.09
## b0 2.17
## b1 0.11
## m0[1] 2.22
## m0[2] 2.69
## m0[3] 2.22
## m0[4] 2.25
## m0[5] 3.00
## m0[6] 3.30
## m0[7] 2.88
## m0[8] 2.68
## m0[9] 3.14
## m0[10] 3.00
## m0[11] 2.88
## m0[12] 3.03
## m0[13] 3.24
## m0[14] 2.39
## m0[15] 2.37
## m0[16] 3.20
## m0[17] 2.77
## m0[18] 3.30
## m0[19] 2.62
## m0[20] 2.51
## y0[1] 9.00
## y0[2] 13.00
## y0[3] 14.00
## y0[4] 5.00
## y0[5] 18.00
## y0[6] 30.00
## y0[7] 20.00
## y0[8] 16.00
## y0[9] 26.00
## y0[10] 21.00
## y0[11] 21.00
## y0[12] 24.00
## y0[13] 30.00
## y0[14] 16.00
## y0[15] 8.00
## y0[16] 27.00
## y0[17] 15.00
## y0[18] 21.00
## y0[19] 19.00
## y0[20] 5.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 657.13 657.44 1.00 0.67 655.08 658.04 1.01 652 574
## b0 2.17 2.17 0.12 0.12 1.97 2.37 1.01 494 394
## b1 0.11 0.11 0.02 0.02 0.09 0.14 1.01 604 619
## m0[1] 2.22 2.22 0.12 0.12 2.03 2.41 1.01 494 388
## m0[2] 2.69 2.69 0.06 0.06 2.59 2.79 1.01 638 611
## m0[3] 2.22 2.23 0.12 0.12 2.03 2.41 1.01 494 387
## m0[4] 2.25 2.25 0.11 0.11 2.06 2.43 1.01 493 388
## m0[5] 2.99 3.00 0.06 0.05 2.90 3.09 1.00 1886 1332
## m0[6] 3.30 3.30 0.08 0.08 3.16 3.43 1.00 1371 1383
## m0[7] 2.88 2.88 0.05 0.05 2.79 2.97 1.00 1283 1363
## m0[8] 2.67 2.67 0.06 0.06 2.57 2.78 1.01 616 512
## m0[9] 3.13 3.14 0.07 0.06 3.02 3.24 1.00 1828 1466
## m0[10] 3.00 3.00 0.06 0.05 2.90 3.09 1.00 1897 1332
## m0[11] 2.88 2.88 0.05 0.05 2.79 2.97 1.00 1293 1363
## m0[12] 3.03 3.03 0.06 0.06 2.93 3.12 1.00 1972 1633
## m0[13] 3.24 3.24 0.08 0.07 3.11 3.36 1.00 1516 1424
## m0[14] 2.39 2.39 0.09 0.09 2.24 2.55 1.01 495 365
## m0[15] 2.37 2.37 0.10 0.09 2.21 2.53 1.01 494 370
## m0[16] 3.19 3.20 0.07 0.07 3.07 3.30 1.00 1657 1445
## m0[17] 2.77 2.77 0.06 0.06 2.67 2.86 1.01 786 764
## m0[18] 3.30 3.30 0.08 0.08 3.16 3.43 1.00 1365 1383
## m0[19] 2.62 2.62 0.07 0.07 2.50 2.73 1.01 561 483
## m0[20] 2.51 2.51 0.08 0.08 2.38 2.64 1.01 511 393
## y0[1] 9.19 9.00 3.15 2.97 4.00 15.00 1.00 1746 1725
## y0[2] 14.73 15.00 3.94 4.45 9.00 21.00 1.00 1923 2035
## y0[3] 9.26 9.00 3.23 2.97 4.00 15.00 1.00 1612 1931
## y0[4] 9.54 9.00 3.30 2.97 5.00 16.00 1.00 1426 1499
## y0[5] 19.92 20.00 4.69 4.45 13.00 28.00 1.00 1929 1807
## y0[6] 27.22 27.00 5.64 5.93 18.00 37.00 1.00 2115 2037
## y0[7] 17.80 18.00 4.38 4.45 11.00 25.00 1.00 1948 2025
## y0[8] 14.57 14.00 4.00 4.45 8.00 22.00 1.00 1825 2074
## y0[9] 22.90 23.00 4.87 4.45 15.00 31.00 1.00 1764 1969
## y0[10] 20.00 20.00 4.54 4.45 13.00 28.00 1.00 1969 1879
## y0[11] 17.85 18.00 4.34 4.45 11.00 25.00 1.00 1888 1972
## y0[12] 20.76 20.00 4.75 4.45 13.00 29.00 1.00 2076 1848
## y0[13] 25.79 26.00 5.35 5.93 17.00 35.00 1.00 1922 1837
## y0[14] 11.03 11.00 3.45 2.97 6.00 17.00 1.00 1682 1927
## y0[15] 10.62 10.00 3.31 2.97 6.00 16.00 1.00 1604 1925
## y0[16] 24.42 24.00 5.38 5.93 16.00 33.00 1.00 1916 1892
## y0[17] 15.96 16.00 4.18 4.45 9.95 23.00 1.00 1915 2002
## y0[18] 27.17 27.00 5.91 5.93 18.00 37.00 1.00 1537 1929
## y0[19] 13.73 14.00 3.76 4.45 8.00 20.00 1.00 1556 1721
## y0[20] 12.49 12.00 3.65 2.97 7.00 19.00 1.00 1547 1854
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
ex13-1.stan
data{
int N;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
real<lower=0> sr0;
vector[N] r0;
}
model{
r0~normal(0,sr0);
for(i in 1:N)
y~poisson_log(b0+r0[i]+b1*x);
}
generated quantities{
vector[N] m0;
array[N] int y0;
for(i in 1:N){
m0[i]=b0+r0[i]+b1*x[i];
y0[i]=poisson_log_rng(m0[i]);
}
}mdl=cmdstan_model('./ex13-1.stan')
data=list(N=n,x=x,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -77904
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 13247.1 0.0196361 261.206 0.4196 0.4196 141
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 148 13353.4 0.00102166 6298 0.5652 0.5652 221
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 13353.40
## b0 1.73
## b1 0.17
## sr0 0.00
## r0[1] 0.00
## r0[2] 0.00
## r0[3] 0.00
## r0[4] 0.00
## r0[5] 0.00
## r0[6] 0.00
## r0[7] 0.00
## r0[8] 0.00
## r0[9] 0.00
## r0[10] 0.00
## r0[11] 0.00
## r0[12] 0.00
## r0[13] 0.00
## r0[14] 0.00
## r0[15] 0.00
## r0[16] 0.00
## r0[17] 0.00
## r0[18] 0.00
## r0[19] 0.00
## r0[20] 0.00
## m0[1] 1.81
## m0[2] 2.52
## m0[3] 1.81
## m0[4] 1.84
## m0[5] 2.98
## m0[6] 3.45
## m0[7] 2.81
## m0[8] 2.50
## m0[9] 3.20
## m0[10] 2.99
## m0[11] 2.81
## m0[12] 3.04
## m0[13] 3.36
## m0[14] 2.07
## m0[15] 2.03
## m0[16] 3.29
## m0[17] 2.64
## m0[18] 3.45
## m0[19] 2.41
## m0[20] 2.24
## y0[1] 9.00
## y0[2] 15.00
## y0[3] 4.00
## y0[4] 6.00
## y0[5] 19.00
## y0[6] 26.00
## y0[7] 11.00
## y0[8] 12.00
## y0[9] 24.00
## y0[10] 20.00
## y0[11] 8.00
## y0[12] 20.00
## y0[13] 30.00
## y0[14] 5.00
## y0[15] 7.00
## y0[16] 32.00
## y0[17] 8.00
## y0[18] 39.00
## y0[19] 13.00
## y0[20] 15.00mcmc=goMCMC(mdl,data,wrm=500)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1 finished in 0.9 seconds.
## Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 2 finished in 0.9 seconds.
## Chain 3 finished in 0.9 seconds.
## Chain 4 finished in 0.9 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.9 seconds.
## Total execution time: 1.0 seconds.seeMCMC(mcmc,'[m,r]')## variable mean median sd mad q5 q95 rhat ess_bulk
## lp__ 13242.54 13242.20 14.49 15.57 13220.30 13268.10 1.05 53
## b0 2.17 2.17 0.03 0.03 2.12 2.22 1.01 652
## b1 0.11 0.11 0.00 0.00 0.11 0.12 1.01 671
## sr0 0.01 0.01 0.01 0.01 0.00 0.02 1.05 55
## r0[1] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 2669
## r0[2] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 3337
## r0[3] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 2843
## r0[4] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3511
## r0[5] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 3200
## r0[6] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 2753
## r0[7] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3004
## r0[8] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3743
## r0[9] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3754
## r0[10] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 4534
## r0[11] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 2822
## r0[12] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 2662
## r0[13] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 2919
## r0[14] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 3858
## r0[15] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3170
## r0[16] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 3022
## r0[17] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3075
## r0[18] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3845
## r0[19] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3313
## r0[20] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3432
## m0[1] 2.22 2.22 0.03 0.03 2.18 2.27 1.01 815
## m0[2] 2.69 2.69 0.02 0.02 2.66 2.72 1.00 1350
## m0[3] 2.23 2.23 0.03 0.03 2.18 2.27 1.00 810
## m0[4] 2.25 2.25 0.03 0.03 2.20 2.29 1.01 797
## m0[5] 3.00 3.00 0.02 0.01 2.97 3.02 1.00 2272
## m0[6] 3.30 3.30 0.02 0.02 3.26 3.34 1.00 1098
## m0[7] 2.88 2.88 0.02 0.01 2.86 2.91 1.00 2001
## m0[8] 2.68 2.68 0.02 0.02 2.65 2.71 1.00 1267
## m0[9] 3.14 3.14 0.02 0.02 3.11 3.17 1.00 2138
## m0[10] 3.00 3.00 0.02 0.02 2.97 3.03 1.00 2745
## m0[11] 2.88 2.88 0.02 0.02 2.86 2.91 1.00 2114
## m0[12] 3.03 3.03 0.02 0.02 3.01 3.06 1.00 2245
## m0[13] 3.24 3.24 0.02 0.02 3.21 3.28 1.00 1705
## m0[14] 2.39 2.39 0.02 0.03 2.35 2.43 1.01 810
## m0[15] 2.37 2.37 0.02 0.02 2.33 2.41 1.01 847
## m0[16] 3.20 3.20 0.02 0.02 3.16 3.23 1.00 1644
## m0[17] 2.77 2.77 0.02 0.02 2.74 2.80 1.00 1497
## m0[18] 3.30 3.30 0.02 0.02 3.27 3.34 1.00 1634
## m0[19] 2.62 2.62 0.02 0.02 2.59 2.65 1.00 1086
## m0[20] 2.51 2.51 0.02 0.02 2.48 2.54 1.00 906
## y0[1] 9.28 9.00 2.99 2.97 5.00 14.00 1.00 2112
## y0[2] 14.73 15.00 3.92 4.45 9.00 22.00 1.00 2117
## y0[3] 9.23 9.00 3.09 2.97 5.00 15.00 1.00 1902
## y0[4] 9.49 9.00 3.01 2.97 5.00 15.00 1.00 1986
## y0[5] 20.06 20.00 4.52 4.45 13.00 28.00 1.00 2115
## y0[6] 27.27 27.00 5.23 5.93 19.00 36.00 1.00 1904
## y0[7] 17.82 18.00 4.17 4.45 11.00 25.00 1.00 1979
## y0[8] 14.57 15.00 3.78 4.45 9.00 21.00 1.00 2096
## y0[9] 22.90 23.00 4.89 4.45 15.00 31.00 1.00 1798
## y0[10] 20.20 20.00 4.56 4.45 13.00 28.00 1.00 1876
## y0[11] 17.87 18.00 4.30 4.45 11.00 25.00 1.00 1946
## y0[12] 20.88 21.00 4.65 4.45 13.00 29.00 1.00 1852
## y0[13] 25.59 25.00 5.25 4.45 17.00 34.00 1.00 1675
## y0[14] 10.98 11.00 3.29 2.97 6.00 17.00 1.00 1874
## y0[15] 10.73 11.00 3.30 2.97 6.00 16.00 1.00 2061
## y0[16] 24.45 24.00 5.12 4.45 16.00 33.00 1.00 2019
## y0[17] 15.77 16.00 3.96 4.45 9.00 23.00 1.00 2034
## y0[18] 27.10 27.00 5.30 5.93 19.00 36.00 1.00 1937
## y0[19] 13.75 13.00 3.78 4.45 8.00 20.00 1.00 1851
## y0[20] 12.36 12.00 3.51 2.97 7.00 18.00 1.00 1943
## ess_tail
## 87
## 1295
## 1172
## 72
## 901
## 741
## 808
## 793
## 788
## 599
## 839
## 757
## 827
## 829
## 674
## 819
## 763
## 647
## 646
## 719
## 649
## 915
## 762
## 736
## 1418
## 1410
## 1269
## 1324
## 1108
## 1333
## 1304
## 1249
## 1199
## 1345
## 1304
## 1299
## 1336
## 1236
## 1170
## 1305
## 975
## 1364
## 1320
## 1329
## 1918
## 2011
## 1853
## 2089
## 1829
## 1755
## 1758
## 1924
## 1649
## 1948
## 1917
## 1983
## 1854
## 1910
## 1972
## 2051
## 1877
## 1711
## 2001
## 1933
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')